homework 6, version 0

Submission by: Jazzy Doe (jazz@mit.edu)

Homework 6: Epidemic modeling III

18.S191, fall 2020

This notebook contains built-in, live answer checks! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.

For MIT students: there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.

Feel free to ask questions!

Let's create a package environment:

We began this module with data on the COVID-19 epidemic, but then looked at mathematical models. How can we make the connection between data and models?

Models have parameters, such as the rate of recovery from infection. Where do the parameter values come from? Ideally we would like to extract them from data. The goal of this homework is to do this by fitting a model to data.

For simplicity we will use data that generated from the spatial model in Homework 5, rather than real-world data, and we will fit the simplest SIR model. But the same ideas apply more generally.

There are many ways to fit a function to data, but all must involve some form of optimization, usually minimization of a particular function, a loss function; this is the basis of the vast field of machine learning.

The loss function is a function of the model parameters; it measures how far the model output is from the data, for the given values of the parameters.

We emphasise that this material is pedagogical; there is no suggestion that these specific techniques should be used actual calculations; rather, it is the underlying ideas that are important.

Exercise 1: Calculus without calculus

Before we jump in to simulating the SIR equations, let's experiment with a simple 1D function. In calculus, we learn techniques for differentiating and integrating symbolic equations, e.g. $\frac{d}{dx}{x}^n=n{x}^{n-1}$. But in real applications, it is often impossible to apply these techniques, either because the problem is too complicated to solve symbolically, or because our problem has no symbolic expression, like when working with experimental results.

Instead, we use ✨ computers ✨ to approximate derivatives and integrals. Instead of applying rules to symbolic expressions, we use much simpler strategies that only use the output values of our function.

As a first example, we will approximate the derivative of a function. Our method is inspired by the analytical definition of the derivative!

$$f^{\prime }(a):=\underset{h\to 0}{lim}\frac{f(a+h)-f(a)}{h}.$$

The finite difference method simply fixes a small value for $h$, say $h={10}^{-3}$, and then approximates the derivative as:

$$f^{\prime }(a)\simeq \frac{f(a+h)-f(a)}{h}.$$

Exercise 1.1 - tangent line

👉 Write a function finite_difference_slope that takes a function f and numbers a and h. It returns the slope $f^{\prime }(a)$, approximated using the finite difference formula above.

👉 Write a function tangent_line that takes the same arguments f, a and g, but it returns a function. This function ($\mathbb{R}\to \mathbb{R}$) is the tangent line with slope $f^{\prime }(a)$ (computed using finite_difference_slope) that passes through $(a,f(a))$.

The slider below controls $h$ using a log scale. In the (mathematical) definition of the derivative, we take $\underset{h\to 0}{lim}$. This corresponds to moving the slider to the left.

Notice that, as you decrease $h$, the tangent line gets more accurate, but what happens if you make $h$ too small?

Exercise 1.2 - antiderivative

In the finite differences method, we approximated the derivative of a function:

$$f^{\prime }(a)\simeq \frac{f(a+h)-f(a)}{h}$$

We can do something very similar to approximate the 'antiderivate' of a function. Finding the antiderivative means that we use the slope $f^{\prime }$ to compute $f$ numerically!

This antiderivative problem is illustrated below. The only information that we have is the slope at any point $a\in \mathbb{R}$, and we have one initial value, $f(1)$.

Using only this information, we want to reconstruct $f$.

By rearranging the equation above, we get the Euler method:

$$f(a+h)\simeq h{f}^{\prime }(a)+f(a)$$

Using this formula, we only need to know the value $f(a)$ and the slope $f^{\prime }(a)$ of a function at $a$ to get the value at $a+h$. Doing this repeatedly can give us the value at $a+2h$, at $a+3h$, etc., all from one initial value $f(a)$.

👉 Write a function euler_integrate_step that applies this formula to a known function $f^{\prime }$ at $a$, with step size $h$ and the initial value $f(a)$. It returns the next value, $f(a+h)$.

👉 Write a function euler_integrate that takes takes a known function $f^{\prime }$, the initial value $f(a)$ and a range T with a == first(T) and h == step(T). It applies the function euler_integrate_step repeatedly, once per entry in T, to produce the sequence of values $f(a+h)$, $f(a+2h)$, etc.

Let's try it out on $f^{\prime }(x)=3x^2$ and T ranging from $0$ to $10$.

We already know the analytical solution $f(x)=x^3$, so the result should be an array going from (approximately) 0.0 to 1000.0.

You see that our numerical antiderivate is not very accurate, but we can get a smaller error by choosing a smaller step size. Try it out!

There are also alternative integration methods that are more accurate with the same step size. Some methods also use the second derivative, other methods use multiple steps at once, etc.! This is the study of Numerical Methods.

Exercise 2: Simulating the SIR differential equations

Recall from the lectures that the ordinary differential equations (ODEs) for the SIR model are as follows:

$$\begin{array}{rl}\dot{s}& =-\beta si\\ \dot{i}& =+\beta si-\gamma i\\ \dot{r}& =+\gamma i\end{array}$$

where $\dot{s}:=\frac{ds}{dt}$ is the derivative of $s$ with respect to time. Recall that $s$ denotes the proportion (fraction) of the population that is susceptible, a number between $0$ and $1$.

We will use the simplest possible method to simulate these, namely the Euler method. The Euler method is not always a good method to solve ODEs accurately, but for our purposes it is good enough.

In the previous exercise, we introduced the euler method for a 1D function, which you can see as an ODE that only depends on time. For the SIR equations, we have an ODE that only depends on the previous value, not on time, and we have 3 equations instead of 1.

The solution is quite simple, we apply the euler method to each of the differential equations within a single time step to get new values for each of $s$, $i$ and $r$ at the end of the time step in terms of the values at the start of the time step. The euler discretised equations are:

$$\begin{array}{rl}s(t+h)& =s(t)-h\cdot \beta s(t)i(t)\\ i(t+h)& =i(t)+h\cdot (\beta s(t)i(t)-\gamma i(t))\\ r(t+h)& =r(t)+h\cdot \gamma i(t)\end{array}$$

👉 Implement a function euler_SIR_step(β, γ, sir_0, h) that performs a single Euler step for these equations with the given parameter values and initial values, with a step size $h$.

sir_0 is a 3-element vector, and you should return a new 3-element vector with the values after the timestep.

👉 Implement a function euler_SIR(β, γ, sir_0, T) that applies the previously defined function over a time range $T$.

You should return a vector of vectors: a 3-element vector for each point in time.

Let's plot $s$, $i$ and $r$ as a function of time.

👉 Do you see an epidemic outbreak (i.e. a rapid growth in number of infected individuals, followed by a decline)? What happens after a long time? Does everybody get infected?

👉 Make an interactive visualization in which you vary $\beta$ and $\gamma$. What relation should $\beta$ and $\gamma$ have for an epidemic outbreak to occur?

Exercise 3: Numerical gradient

For fitting we need optimization, and for optimization we will use derivatives (rates of change). In Exercise 1, we wrote a function finite_difference_slope(f, a) to approximate $f^{\prime }(a)$. In this exercise we will write a function to compute partial derivatives.

Exercise 3.1

👉 Write functions ∂x(f, a, b) and ∂y(f, a, b) that calculate the partial derivatives $\frac{\partial f}{\partial x}$ and $\frac{\partial f}{\partial y}$ at $(a,b)$ of a function $f:{\mathbb{R}}^2\to \mathbb{R}$ (i.e. a function that takes two real numbers and returns one real).

Recall that $\frac{\partial f}{\partial x}$ is the derivative of the single-variable function $g(x):=f(x,b)$ obtained by fixing the value of $y$ to $b$.

You should use anonymous functions for this. These have the form x -> x^2, meaning "the function that sends $x$ to $x^2$".

Exercise 3.2

👉 Write a function gradient(f, a, b) that calculates the gradient of a function $f$ at the point $(a,b)$, given by the vector $\mathrm{\nabla }f(a,b):=(\frac{\partial f}{\partial x}(a,b),\frac{\partial f}{\partial y}(a,b))$.

Exercise 4: Minimisation using gradient descent

In this exercise we will use gradient descent to find local minima of (smooth enough) functions.

To do so we will think of a function as a hill. To find a minimum we should "roll down the hill".

Exercise 4.1

We want to minimize a 1D function, i.e. a function $f:\mathbb{R}\to \mathbb{R}$. To do so we notice that the derivative tells us the direction in which the function increases. Positive slope means that the minimum is to the left, negative slope means to the right. So our gradient descent method is to take steps in the opposite direction, of a small size $\eta \cdot f^{\prime }(x_0)$.

👉 Write a function gradient_descent_1d_step(f, x0) that performs a single gradient descrent step, from the point x0 and using your function finite_difference_slope to approximate the derivative. The result should be the next guess for $x$.

$x_0=$

-1.0

👉 Write a function gradient_descent_1d(f, x0) that repeatedly applies the previous function (N_steps times), starting from the point x0, like in the vizualisation above. The result should be the final guess for $x$.

Right now we take a fixed number of steps, even if the minimum is found quickly. What would be a better way to decide when to end the function?

Exericse 4.2

Multivariable calculus tells us that the gradient $\mathrm{\nabla }f(a,b)$ at a point $(a,b)$ is the direction in which the function increases the fastest. So again we should take a small step in the opposite direction. Note that the gradient is a vector which tells us which direction to move in the plane $(a,b)$. We multiply this vector with the scalar $\eta$ to control the step size.

👉 Write functions gradient_descent_2d_step(f, x0, y0) and gradient_descent_2d(f, x0, y0) that do the same for functions $f(x,y)$ of two variables.

$x_0=$

0.0

$y_0=$

0.0

We also prepared a 3D visualisation if you like! It's a bit slow...

👉 Can you find different minima?

Exercise 5: Learning parameter values

In this exercise we will apply gradient descent to fit a simple function $y=f_{\alpha ,\beta }(x)$ to some data given as pairs $(x_i,y_i)$. Here $\alpha$ and $\beta$ are parameters that appear in the form of the function $f$. We want to find the parameters that provide the best fit, i.e. the version $f_{\alpha ,\beta }$ of the function that is closest to the data when we vary $\alpha$ and $\beta$.

To do so we need to define what "best" means. We will define a measure of the distance between the function and the data, given by a loss function, which itself depends on the values of $\alpha$ and $\beta$. Then we will minimize the loss function over $\alpha$ and $\beta$ to find those values that minimize this distance, and hence are "best" in this precise sense.

The iterative procedure by which we gradually adjust the parameter values to improve the loss function is often called machine learning or just learning, since the computer is "discovering" information in a gradual way, which is supposed to remind us of how humans learn. [Hint: This is not how humans learn.]

Exercise 5.1 - 🎲 frequencies

We generate a small dataset by throwing 10 dice, and counting the sum. We repeat this experiment many times, giving us a frequency distribution in a familiar shape.

Let's try to fit a gaussian (normal) distribution. Its PDF with mean $\mu$ and standard deviation $\sigma$ is

$$f_{\mu ,\sigma }(x):=\frac{1}{\sigma \sqrt{2\pi }}\exp [-\frac{(x-\mu {)}^2}{2{\sigma }^2}]$$

👉 (Not graded) Manually fit a Gaussian distribution to our data by adjusting $\mu$ and $\sigma$ until you find a good fit.

μ =

24.0

σ =

12.0

Show manual fit:

What we just did was adjusting the function parameters until we found the best possible fit. Let's automate this process! To do so, we need to quantify how good or bad a fit is.

👉 Define a loss function to measure the "distance" between the actual data and the function. It will depend on the values of $\mu$ and $\sigma$ that you choose:

$$\mathcal{L}(\mu ,\sigma ):=\sum _i[f_{\mu ,\sigma }(x_i)-y_i{]}^2$$

👉 Use your gradient_descent_2d function to find a local minimum of $\mathcal{L}$, starting with initial values $\mu =30$ and $\sigma =1$. Call the found parameters found_μ and found_σ.

Go back to the graph to see your optimized gaussian curve!

If your fit is close, then probability theory tells us that the found parameter $\mu$ should be close to the weighted mean of our data, and $\sigma$ should approximate the sample standard deviation. We have already computed these values, and we check how close they are:

Exercise 6: Putting it all together — fitting an SIR model to data

In this exercise we will fit the (non-spatial) SIR ODE model from Exercise 1 to some data generated from the spatial model in Problem Set 4. If we are able to find a good fit, that would suggest that the spatial aspect "does not matter" too much for the dynamics of these models. If the fit is not so good, perhaps there is an important effect of space. (As usual in statistics, and indeed in modelling in general, we should be very cautious of making claims of this nature.)

This fitting procedure will be different from that in Exercise 4, however: we no longer have an explicit form for the function that we are fitting – rather, it is the output of an ODE! So what should we do?

We will try to find the parameters $\beta$ and $\gamma$ for which the output of the ODEs when we simulate it with those parameters best matches the data!

Exercise 6.1

Below the result from Homework 4, Exercise 3.2. These are the average S, I, R fractions of running 20 simulations. Click on it!

👉 (Not graded) Manually fit the SIR curves to our data by adjusting $\beta$ and $\gamma$ until you find a good fit.

β =

0.05

γ =

0.005

Show manual fit:

👉 To do this automatically, we will again need to define a loss function $\mathcal{L}(\beta ,\gamma )$. This will compare the solution of the SIR equations with parameters $\beta$ and $\gamma$ with your data.

This time, instead of comparing two vectors of numbers, we need to compare two vectors of vectors, the S, I, R values.

👉 Use this loss function to find the optimal parameters $\beta$ and $\gamma$.

Function library

Just some helper functions used in the notebook.